C............................<HMF2.FOR>........... 5-MAY-1986 11:41
      SUBROUTINE HMF2(I,Z,ZKM,OX,N2,O2,TN,TE,UM,DIP,CHI,OXJ,OPLS
     > ,OSAV,DT)
C....... This subroutine uses a Newton procedure to solve the equation
C....... Beta=D*sinI*sinI/H/H-U*CosI*sinI/H for HmF2
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      IF(DIP.EQ.0) DIP=0.1
      U=UM*100
      Z0=ZKM*1.0E5
      Z=Z0
      JTER=0
      H=0.0001*Z
C
C....... start Newton
 20   CALL EFFZ(I,Z,Z0,OX,N2,O2,TN,TE,U,DIP,FZ)
      FEX=FZ
      Z=Z+H
      CALL EFFZ(I,Z,Z0,OX,N2,O2,TN,TE,U,DIP,FZ)
      DEX=(FZ-FEX)/H
      IF(DEX.EQ.0) GO TO 30
      Z=Z-H-FEX/DEX
      JTER=JTER+1
      IF(JTER.GT.22) GO TO 30
      IF(ABS(FEX/(DEX*Z)).GT.1.0E-3) GO TO 20
      CALL OPLSD(I,Z,Z0,OX,N2,O2,TN,TE,OXJ,OPLS,OSAV,DT,CHI,DIP)
      Z=Z*1.0E-5
      RETURN
 30   CONTINUE
      PRINT*, '  *** JTER > 22'
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE EFFZ(I,Z,Z0,OX,N2,O2,TN,TE,U,DIP,FZ)
C....... This subroutine Contains the funCtion to be solved by Newton
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
C....... evaluate 1/H where H=sCale height, and density at altitude z
      CALL NATM(Z,Z0,OX,N2,O2,TN,RHOX,RHN2,RHO2,OXZ,N2Z,O2Z)
C....... evaluate reaCtion rates
      CALL RRATES(I,TN,TN,TE,RN2,RO2)
C....... evaluate [O+] loss rates
      BETA=RN2*N2Z+RO2*O2Z
C....... evaluate O+ - O Collision frequenCy
      NU=6.82E-10*N2Z+3.67E-11*OXZ*SQRT(TN)*(1-0.064*DLOG10(TN))**2
C....... diffusion CoeffiCient
      DH=RHOX
      DIFF=5.15E6*(TE+TN)/NU
      FZ=BETA-DIFF*(SIN(DIP)*DH)**2-U*SIN(DIP)*COS(DIP)*DH
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE NATM(Z,Z0,OX,N2,O2,TN,RHOX,RHN2,RHO2,OXZ,N2Z,O2Z)
C........ evaluate reCiproCal of sCale heights and densities at z
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      GRAV=3.98E+20/(6.371E+8+Z)**2
      RHOX=1.94E-7*GRAV/TN
      RHN2=1.75*RHOX
      RHO2=2.0*RHOX
      OXZ=OX*EXP(-(Z-Z0)*RHOX)
      N2Z=N2*EXP(-(Z-Z0)*RHN2)
      O2Z=O2*EXP(-(Z-Z0)*RHO2)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE RRATES(I,TN,TI,TE,RN2,RO2)
C....... [O+] reaCtion rates
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
C......... O+ + N2 --> NO+ + N St. MauriCe and Torr JGR ,1978, p969
      T13=(4.*TN+7.*TI)/11.0/300
      IF(T13.LE.5.67) RN2=1.533E-12-5.92E-13*T13+8.60E-14*T13**2
      IF(T13.GT.5.67) RN2=2.73E-12-1.15E-12*T13+1.483E-13*T13**2
C
C........ O+ + O2 --> O2+ + 0 .  from Chen et al j.Chem.phys. 1979
C     RO2=2.1E-11*((TN+2*TI)/3/300)** -0.763
      RO2=2.1E-11*((TN+2*TI)/900)**(-0.763)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE OPLSD(I,Z,Z0,OX,N2,O2,TN,TE,OXJ,OPLS,OSAV,DT,CHI,DIP)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      DIMENSION XN(3),COLUMN(3),SIGABS(3)
      DATA SIGABS/8E-18,11E-18,11E-18/
C
C......... Evaluate [O+] at peak
      CALL NATM(Z,Z0,OX,N2,O2,TN,RHOX,RHN2,RHO2,OXZ,N2Z,O2Z)
      XN(1)=OXZ
      XN(2)=O2Z
      XN(3)=N2Z
      CALL COLUM(IJ,CHI,Z,TN,XN,COLUMN)
      TAU=0.0
      DO 4 IS=1,3
 4    TAU=TAU+SIGABS(IS)*COLUMN(IS)
C
      CALL RRATES(I,TN,TN,TE,RN2,RO2)
      BETA=RN2*N2Z+RO2*O2Z
C...... PROD/LOSS. 1+1.75*SIN(DIP) produCtion due to diffusion
      POBETA=(1+1.75*SIN(DIP)) *OXJ*OXZ*EXP(-TAU)/BETA
      OPLS=POBETA+(OSAV-POBETA)*DEXP(-BETA*DT)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE COLUM(J,CHI,Z,TN,XN,COLUMN)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C++++ this routine evaluates the neutral Column density
C++++ see smith & smith jgr 1972 p 3592, subr ambs is the msis
C++++ neutral atmosphere model for z<120 km, Chi=solar zenith
C++++ angle, re & ge radius and grav Con for earth
      DIMENSION XN(3),COLUMN(3),SN(5),M(3)
      DATA A,B,C,D,F,G/1.0606963,0.55643831,1.0619896,1.7245609
     >  ,0.56498823,0.06651874/
        DATA    EM   ,   M(1) , M(2) , M(3) ,  RE   , GE
     1    / 1.662E-24 ,   16. ,  32. ,  28. ,6.357E8, 980/
C------ is sza>90.0 degrees
      IF(CHI.LT.1.5708)GO  TO 2938
C------ CalCulate grazing inCidenCe parameters
      ALTG=(6371.0E5+Z)*SIN(3.1416-CHI)-6371.0E5
      IF(ALTG.GE.85.E5)GO TO 1001
C------ if altg < 85km, produCtion=0, make Column large to aChieve this
      DO 20 I=1,3
 20   COLUMN(I)=1.E+25
      RETURN
C
C.... Call neutral atmosphere for density at grazing inCidenCe
 1001  ZG=ALTG*1.E-5
       RHOX=1.94E-7* 3.98E+20/(6.371E+8+Z)**2/TN
       SN(1)=XN(1)*EXP(-(ALTG-Z)*RHOX)
       SN(2)=XN(2)*EXP(-(ALTG-Z)*RHOX*2.0)
       SN(3)=XN(3)*EXP(-(ALTG-Z)*RHOX*1.75)
C---- to pt p, sh=sCale height, rg=distanCe to pt g, hg=sCale height at g
C
 2938 CONTINUE
      GR=GE*(RE/(RE+Z))**2
      RP=RE+Z
      DO 10 I=1,3
      SH=(1.38E-16*TN)/(EM*M(I)*GR)
      XP=RP/SH
      Y=SQRT(0.5*XP)*ABS(COS(CHI))
      IF(Y.GT.100.0) WRITE(3,100)
  100 FORMAT('WARNING, Y IN COLUMN(I) > 100')
      IF(Y.GT.8) ERFY2=F/(G+Y)
      IF(Y.LT.8) ERFY2=(A+B*Y)/(C+D*Y+Y*Y)
    4 IF(CHI.GT.1.5708)GO  TO 2
      CHAPFN=SQRT(0.5*3.1416*XP)*ERFY2
      COLUMN(I)=XN(I)*SH*CHAPFN
        GO TO 10
    2 RG=RP*SIN(3.1416-CHI)
      HG=1.38E-16*TN/
     1    (EM*M(I)*980.*(6371.E5/(6371.E5+ALTG))**2)
      XG=RG/HG
      COLUMN(I)=SQRT(0.5*3.1416*XG)*HG*(2.0*SN(I)-XN(I)*ERFY2)
10       CONTINUE
      RETURN
      END